home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 007 / erf.arc / ERF.PLB next >
Encoding:
Text File  |  1985-08-12  |  6.5 KB  |  160 lines
  1. {     ERF.PLB #2.00 85-08-06 GAUSSIAN ERROR-FUNCTION APPROXIMATION
  2.  
  3.        V02 L00 85-08-06 conversion to Turbo Pascal by Dennis E. Hamilton,
  4.                for separate testing for use as part of the Chi-Squared
  5.                Distribution procedure, Q2(x, df).
  6.  
  7.        V01 L01 79-03-22 implementated in Benton Harbor BASIC 05.01.01 by DEH;
  8.                used in testing that BASIC's awful random-number generator.
  9.            L00 79-03-18 Texas Instruments SR52 calculator program derived
  10.                by Dennis E. Hamilton to verify basic recurrence method by
  11.                checking against values in standard tables.       }
  12.  
  13.  
  14. function
  15.  
  16.    erf(z: real {limit of integration} )
  17.  
  18.        :real {value of the Gaussian Error Function at z};
  19.  
  20.  
  21.       {The Gaussian error function is the integral, from 0 to z, of
  22.  
  23.                 2/sqrt(pi)*exp(-sqr(t)) dt
  24.  
  25.        with negative z also defined by erf(-z) = -erf(z).
  26.  
  27.        This function can be used as the basis for computation of the related
  28.        Gaussian probability functions, including the Normal Probability
  29.        distribution and the chi-Squared distribution functions. }
  30.  
  31.  
  32.    const
  33.  
  34.       e1max = 1E-10;
  35.               {largest positive value of x for which exp(x) is exactly 1.0
  36.                to the floating-point precision of this computer.  This value
  37.                is used to avoid computation of insignificant terms.}
  38.  
  39.       e0min = 81;
  40.               {smallest nonzero value of x at which exp(-x) underflows to
  41.                0.0 within the computer's floating-point precision.  This value
  42.                is used to steer around avoidable arithmetic overflows.}
  43.  
  44.  
  45.    function H1(x: real)
  46.  
  47.                :real {approximation of 1-erf(abs(x)) where erf(x) is
  48.                       the standard Gaussian error integral [HMF: 7.1.1]};
  49.  
  50.        const   a6 = 0.0000430638;      a5 = 0.0002765672;
  51.                a4 = 0.0001520143;      a3 = 0.0092705272;
  52.                a2 = 0.0422820123;      a1 = 0.0705230784;
  53.  
  54.              xmax = 9 {value of sqrt(e0min)};
  55.  
  56.        begin {Computation to over 6 decimal-place precision based on Hastings
  57.               approximation 7.1.28 in [Abramowitz, Milton., Stegun, Irene A.
  58.               (eds.)  HANDBOOK OF MATHEMATICAL FUNCTIONS.  National Bureau of
  59.               Standards Applied Mathematics Series Publication #55.  Dover
  60.               Publications (New York: 1965)].  The 4-place version given
  61.               there [HMF: 7.1.26] is also perfectly usable, and is to be
  62.               recommended when e1max is larger than 1E-6 or so.  The present
  63.               form has been chosen to avoid reliance on exp(x), so that
  64.               availability and reliability of any built-in exp(x) function
  65.               does not have to be presumed in order to obtain erf(x) and its
  66.               relatives, such as the normal probability function P(x). }
  67.  
  68.        x := abs(x);
  69.        if x > xmax
  70.        then H1 := 0
  71.        else begin
  72.             x := (((((a6*x + a5)*x + a4)*x + a3)*x + a2)*x + a1)*x + 1.0;
  73.             H1 := sqr(sqr(sqr(sqr(1.0/x))));
  74.             end;
  75.        end {H1};
  76.  
  77.  
  78.    BEGIN {erf};
  79.  
  80.    if z < 0.0
  81.    then erf := H1(z) - 1.0
  82.    else erf := 1.0 - H1(z);
  83.  
  84.       {Note that the loss of precision invariably accompanying addition
  85.        and subtraction of 1.0 when H1(z) is small can often be forestalled by
  86.        using H1 directly in places where erf(z) is required and the sign of
  87.        z is already known. }
  88.  
  89.    END {erf};
  90.  
  91. {  The following results were obtained with a test program (ERFT1.PAS)
  92.    compiled and executed using Borland International Turbo Pascal 3.00A for
  93.    CP/M-80.  This version maintains a 40-bit floating-point significand and
  94.    decimal output conversion rounded to 11 significant digits.   The error
  95.    terms represent the absolute difference found by consulting known 
  96.    tabulations provided along with the z values in file ERFT1.DAT:
  97.  
  98.  
  99. ERFT1> #2.00 85-08-06 STANDARD ERROR FUNCTION TEST & DEMONSTRATION
  100.  
  101.             z                erf(z)           abs(error)
  102.  
  103.        0.000000000E+00   0.0000000000E+00      0.00E+00
  104.        2.579182479E-11   0.0000000000E+00      2.58E-11
  105.        2.579282480E-11   2.9103830457E-11      3.31E-12
  106.        1.000000000E-10   8.7311491370E-11      1.27E-11
  107.        0.000001000       0.00000112836         1.28E-07
  108.  
  109.        0.010000000       0.01128332760         8.80E-08
  110.        0.020000000       0.02256441954         1.55E-07
  111.        0.030000000       0.03384101848         2.04E-07
  112.        0.050000000       0.05637172353         2.54E-07
  113.        0.060000000       0.06762133443         2.60E-07
  114.        0.080000000       0.09007788447         2.41E-07
  115.        0.090000000       0.10128037352         2.20E-07
  116.  
  117.        The values computed for erf(z), with z small, have the 
  118.        indicated  variance from table values of  The Handbook
  119.        of Mathematical Functions [HMF: Table 7.1].    Results
  120.        to within 1E-5 are usually quite acceptable.
  121.  
  122.             z                erf(z)           abs(error)
  123.  
  124.        0.180000000       0.20093593110         9.21E-08
  125.        0.370000000       0.39920611672         1.33E-07
  126.  
  127.        0.460000000       0.48465529432         9.57E-08
  128.        0.470000000       0.49374493215         1.19E-07
  129.        0.480000000       0.50274953019         1.41E-07
  130.        0.490000000       0.51166810037         1.61E-07
  131.        0.500000000       0.52049969828         1.80E-07
  132.  
  133.        0.600000000       0.60385583273         2.58E-07
  134.        0.810000000       0.74800336344         8.28E-08
  135.  
  136.        Midrange values, near .5, reveal the best that can be 
  137.        obtained when either erf(z) or  erfc(z) = 1-erf(c) is
  138.        used.  Test values are from [HMF: Table 7.1].
  139.  
  140.             z                erf(z)           abs(error)
  141.  
  142.        1.830000000       0.99034701019         2.05E-07
  143.        2.506628275       0.99960703858         2.11E-07
  144.        3.069980124       0.99998570912         1.46E-07
  145.        3.544907702       0.99999943889         2.59E-08
  146.        3.963327298       0.99999997623         2.97E-09
  147.        4.341607527       0.99999999890         2.95E-10
  148.  
  149.        4.689472100       0.99999999994         5.46E-11
  150.        5.000000000       1.00000000000         1.82E-12
  151.        1.000000000E+01   1.00000000000         0.00E+00
  152.  
  153.        Near the tail, as erf(z) approaches 1.0, [HMF: Table 7.3]
  154.        values of  erfc(z) = 1 - erf(z)  are used to see how well
  155.        accuracy is maintained despite adjustments by 1.0 in erf.
  156.  
  157. ERFT1> End of test.
  158.                                 }
  159.             (* end of ERF.PLB *)
  160.